Use R, leaflet.js, and dygraphs.js to assess hail risk
Author 1, Affiliation
Author 2, Affiliation
Imagine baseball-sized ice balls falling from the sky.
For some Americans, this is a real weather phenomenon that risks the physical well-being of people and property. But how often does it happen? Where do these events normally happen?
As it turns out, the National Oceanographic and Atmospheric Administration (NOAA) collects massive amounts of data on precipitation using radar stations across the country. Algorithms can be built on top of the data to track ‘hail signatures’ detecting when these events occur and the severity. This severe storm data is captured in NOAA's Severe Weather Data Inventory (SWDI) housed within the National Centers for Environmental Information (NCEI).
Looking at nearly 10 million events over the last 10 years, we can see that most of these events happen in the Midwest.
A closer look at severe events (e.g. hail diameter > 3in), the majority of these events occur in a small snaking section of the country indicating that it's relatively isolated.
[Section to be completed] Get Started with SWDI
SWDI data can be used for everything from agriculture to risk management.
setwd("/Users/jeffreychen/Google Drive/DOC/027-SWDI/cviz_tutorial")
library(sqldf)
library(RColorBrewer)
library(leaflet)
library(googleVis)
library(ggplot2)
library(DT)
library(rgdal)
library(sp)
##Set Parameters
start_date = "2005-01-01"
end_date = "2014-12-31"
range <- as.Date(end_date,"%Y-%m-%d")-as.Date(start_date,"%Y-%m-%d")
series = "nx3hail"
fraction = 0.25
#Load in files
shp <- readOGR("cb_2014_us_county_20m", "cb_2014_us_county_20m", stringsAsFactors=FALSE, encoding="UTF-8")
raw <- read.csv("hail_10yr.csv")
fips <- read.delim("http://www2.census.gov/geo/docs/reference/state.txt",sep="|")
fips <- fips[,c("STATE","STUSAB")]
fips$STATE <- as.numeric(fips$STATE)
fips$STUSAB <- as.character(fips$STUSAB)
#Cutdown
raw <- raw[raw$LON<(-50) & raw$LON>(-140) & raw$LAT > 25,]
raw$LON <- round(raw$LON/fraction)*fraction
raw$LAT <- round(raw$LAT/fraction)*fraction
##De-duplicate by day, latitude and longitude
deduped_day <- sqldf("SELECT DATE, LON, LAT, MAXSIZE
FROM raw
GROUP BY DATE, LON, LAT, MAXSIZE")
deduped_day$lvl_2in<-0
deduped_day$lvl_2in[deduped_day$MAXSIZE>2]<-1
deduped_day$lvl_3in<-0
deduped_day$lvl_3in[deduped_day$MAXSIZE>3]<-1
##Daily gridded frequencies
singles <- sqldf("SELECT LON, LAT,COUNT(date) cnt, SUM(lvl_2in) cnt_2, SUM(lvl_3in) cnt_3
FROM deduped_day
GROUP BY LON, LAT")
for(i in 3:ncol(singles)){
singles[,i]<-singles[,i]/as.numeric(range)
}
write.csv(singles,"singles.csv",row.names=F)
##Set up spatial
points<-singles
coordinates(points)=~LON+LAT
proj4string(points)=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
shp <- spTransform(shp, CRS("+proj=longlat +ellps=WGS84 +datum=WGS84") )
##Spatial join
a<-over(points,shp)
singles <- cbind(singles,a[,c("STATEFP","NAME")])
singles$NAME <- as.character(singles$NAME)
#Add state abbreviation via FIPS
singles$STATEFP <-as.numeric(singles$STATEFP)
singles <- merge(singles,fips,by.x="STATEFP",by.y="STATE",all.x=T, sort=F)
singles$loc <- paste("Nearest County: ",singles$NAME,", ",singles$STUSAB,sep="")
singles$loc[is.na(singles$NAME)]<-""
singles <- singles[!is.na(singles$STUSAB),]
#Subset data into separate frames (needed for layers in leaflet.js)
in3 <- singles[singles$cnt_3>0,]
#Vectorize data for inclusion for popup text
#All points
county_all <- singles$loc
x_all <- singles$LON
y_all <- singles$LAT
pr <- round(100*singles$cnt,3)
#Hail balls > 3in
county_3 <- in3$loc
x3 <- in3$LON
y3 <- in3$LAT
pr3 <- round(100*in3$cnt_3,3)
#Popup
content_all <- paste("<h3>",county_all,"</h3>Grid Point: ",x_all,", ",y_all,"<p>Prob of Any Hail <span style='color:red'><strong>: ",pr,"%</strong></span></p>")
content_3 <- paste0("<h3>",county_3,"</h3>Grid Point: ",x3,", ",y3,"<p>Prob of Hail (> 3in)<span style='color:red'><strong>: ",pr,"%</strong></span></p>")
pal <- colorNumeric(
palette = "Set1",
domain = pr
)
pal2 <- colorNumeric(
palette = "Set1",
domain = pr3
)
leaflet(width="100%") %>%
setView(lat = mean(y_all), lng = mean(x_all),4) %>%
addTiles('http://{s}.basemaps.cartocdn.com/dark_all/{z}/{x}/{y}.png') %>%
addCircleMarkers(data = singles, lat = ~ LAT, lng = ~ LON,radius=(pr/5),
fillOpacity = 0.8,stroke = FALSE,
color = ~pal(pr), popup = content_all) %>%
addLegend("bottomright", pal = pal, values = pr,
title = "Prob(Hail)",labFormat = labelFormat(suffix = "%")
)
leaflet(width="100%") %>%
setView(lat = mean(y_all), lng = mean(x_all),4) %>%
addTiles('http://{s}.basemaps.cartocdn.com/dark_all/{z}/{x}/{y}.png') %>%
addCircleMarkers(data = in3, lat = ~ LAT, lng = ~ LON, radius=2*(pr3),
fillOpacity = 0.8,stroke = FALSE,
color = ~pal2(pr3), popup = content_all) %>%
addLegend("bottomright", pal = pal2, values = pr3,
title = "Prob(Hail diameter > 3in)",labFormat = labelFormat(suffix = "%")
)